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Abstract 

We present a novel approach to finding critical points in cell-wise 
barycentrically or bilinearly interpolated vector fields on surfaces. The 
Poincare index of the critical points is determined by investigating the 
qualitative behavior of 0-level sets of the interpolants of the vector field 
components in parameter space using precomputed combinatorial results, 
thus avoiding the computation of the Jacobian of the vector field at the 
critical points in order to determine its index. The locations of the critical 
points within a cell are determined analytically to achieve accurate results. 
This approach leads to a correct treatment of cases with two first-order 
critical points or one second-order critical point of bilinearly interpolated 
vector fields within one cell, which would be missed by examining the lin- 
earized field only. We show that for the considered interpolation schemes 
determining the index of a critical point can be seen as a coloring problem 
of cell edges. A complete classification of all possible colorings in terms of 
the types and number of critical points yielded by each coloring is given 
using computational group theory. We present an efficient algorithm that 
makes use of these precomputed classifications in order to find and classify 
critical points in a cell-by-cell fashion. Issues of numerical stability, con- 
struction of the topological skeleton, topological simplification, and the 
statistics of the different types of critical points are also discussed. 

Keywords: Vector field topology, interpolation, barycentric interpo- 
lation, linear interpolation, bilinear interpolation, level sets, higher-order 
singularities, computational group theory, colorings. 

1 Introduction 



The visualization of vector field topology is a problem that arises naturally when 
studying the qualitative structure of flows that are tangential to some surface. 
As usual, we use the term surface for a real, smooth 2-manifold (equipped with 
an atlas consisting of charts), see for example |Q'N83j for an introduction to 
Riemannian Geometry. Having its roots in the theory of dynamical systems, 
the topological skeleton of a Hamiltonian flow on a surface with isolated critical 
points consists of these critical points and trajectories (streamlines) of the vector 
field that lie at the boundary of a hyperbolic sector and connect two of the 
critical points. Hclman and Hcsselink |HH90| introduced the concept of the 
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(a) (b) (c) 

Figure 1: Classification of critical points for three intersection cases for the 
0-level sets of the two vector field components /i (cyan) and (orange) of 
a bilincarly interpolated vector field / = (/i,/2) T : (a) no intersection of the 
level sets, (b) touching of the level sets yielding one critical point, (c) double 
intersection of the level sets yielding two critical points (one saddle and one non- 
saddle). The colors of the vector arrows encode the types of characteristic areas 

as defined in Section B~2l (grccn=++. ycllow=H — , bluc= — h, rcd= ), a set 

of streamlines is shown in gray, and critical points lying in the intersection set 
of the two level sets are shown as black dots. Each vertex of a square is marked 
with ++, + — , — +, according to the sign of j\ and at that vertex. 



topology of a planar vector field to the visualization community and proposed 
the following construction scheme: (1) critical points are located, (2) classified, 
and then (3) trajectories along hyperbolic sectors are traced and connected to 
their originating and terminating critical points or boundary points. Step (2) — 
the classification of a critical point — is usually based on the Jacobian of the 
vector field, see |Per91j . Trajectories of step (3) are typically constructed by 
solving an ordinary differential equation for particle tracing. 

Although a vast body of previous work in the field of flow visualization 
focuses on the problem of how to extend the method of Hclman and Hcssclink 
to vector fields on arbitrary surfaces as well as the second and third step of 
the above algorithm, not much attention has been paid to the first step. In 
this paper, we specifically address the identification and classification of critical 
points in parameter space. 

As efficient computer-based visualization algorithms usually work with dis- 
crete parametrized versions of the surfaces involved — examples of popular dis- 
cretization schemes are for example triangulated or quadrangulated versions of 
the surface — we will in this paper not focus on the well-researched field of how 
to parametrize a given surface (see |FH05j for a recent survey) but assume that a 
surface always comes equipped with a globally continuous discrete parametriza- 
tion that allows a cell-wise (local) barycentric or bilinear interpolation scheme 
of a vector field tangential to the surface in parameter space. 

While this task is rather easy for linear vector fields, the problem setting be- 
comes more interesting for bilinearly interpolated fields. Bilinear interpolation 
is ubiquitous in scientific visualization because it is popular for widely used uni- 
form or curvilinear grids representing planar or curved surfaces. Since bilinear 
interpolation is not linear, it can lead to higher-order critical points, which are 
neglected by often-used linearization approaches. 

In this paper, we introduce a new method that locates and classifies all crit- 
ical points within piecewise linearly (barycentrically) or bilinearly interpolated 
two-dimensional grids. Our method determines the index of a critical point 
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without the need to evaluate the Jacobian of the vector field in the critical 
point in order to determine its Poincare index. For the case of bilinearly inter- 
polated vector fields, our method is able to detect higher-order critical points 
and the presence of two first-order critical points within one cell, which, to our 
knowledge, has not been achieved with the common methods }HH91| for the 
bilinear interpolation scheme before. Figure [T] shows a corresponding example 
and illustrates our classification method. Our approach is based on the idea that 
investigating the qualitative behavior of 0-lcvel sets of the components of the 
interpolated vector field provides information needed to compute the Poincare 
index of a critical point. All qualitatively different possibilities of this behavior 
and the types of critical points yielded by each possibility are completely classi- 
fied using the computational group theory tool GAP (Groups, Algorithms, and 
Programming) GAP06 . Sec the enumeration of cases for marching cubes and 
generic substitope algorithms by Banks et al. |BLS 04 for a previous example 
of an application of computational group theory in the field of scientific visu- 
alization. Furthermore, we discuss a cell-based topology simplification method 
as well as the question of numerical stability. Our approach results in an effi- 
cient, accurate, and robust cell-based algorithm for detecting all critical points 
of barycentrically or bilinearly interpolated 2D vector fields. 

The paper is organized as follows. First, we will give a short review of the 
visualization literature dealing with vector field topology. Then, the theoretical 
foundations of vector field topology, namely the theory of the qualitative behav- 
ior of second-order dynamical systems along with such fundamental notions as 
those of critical points, separatrices, and the Poincare index of a critical point 
are reviewed. Following this, we will present our new approach — first the gen- 
eral framework will be discussed and then applied to two cases: barycentrically 
and bilinearly interpolated vector fields. Then, we will deal with open issues 
such as critical points on the boundary of cells and numeric stability followed 
by a more detailed description of the cell-based algorithm. We will conclude 
giving results and a short review of our method. 

This paper has accompanying material in the form of online resources, 
namely the GAP programs used in this paper (Online Resource 3) and lists 
of equivalence classes of colored cells referred to in Theorem 15.31 (Online Re- 
source 1) and Theorem 16.71 (Online Resource 2). 



2 Previous Work 

Topology-based methods for planar vector fields were first proposed to the visu- 
alization community by Helman and Hesselink I II !!)(")■ , employing methods from 
the theory of dynamical systems |ALGM73j to planar linear (or linearized) vec- 
tor fields in order to visualize flow characteristics. Their work has triggered a 
large body of further research on topology-based flow visualization, an overview 
of which is given by Post et al. |PVH+03] and Scheuermann and Tricoche jST05j . 



The case of a planar linear vector field is relatively simple to deal with, but it 
has a couple of drawbacks, most notably that only first-order critical points can 
be detected. Much effort has been put into methods to overcome this drawback 
imposed by the interpolation scheme and to addresses the issue of detecting 
and processing higher-order critical points of interpolated vector fields, also of 
vector fields on arbitrary surfaces jLCJK + 09) . Such higher-order critical points 
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can be found, for example, by using piecewise linear interpolation schemes in 
combination with a clustering of first-order critical points according to Tricoche 
et al. jTSHOO] . Another example is the work by Theisel |The02j , who proposes 
a method for designing piecewise linear planar vector fields of arbitrary topol- 
ogy. Nonlinear interpolation schemes have been investigated by Scheuermann et 
al. [SHK+97] and Zhang et al. |ZMT06j . Scheuermann et al. 1SKMR98) propose 
a way to approximate higher-order critical points using Clifford algebra. Li ct 
al. |LVRL06] use interpolation schemes based on a polar coordinate representa- 
tion to detect vector field singularities on a surface. 

In contrast to the global variational approach taken in [PP03] in which the 
authors construct a discrete Hodge decomposition in order to obtain location 
and type of critical points of vector fields on polyhedral surfaces, our approach 
is local and cell-oriented. It may thus be easier to implement when a discrete 
parametrization of the surface is already given and can be used in conjunction 
with other local, grid-based methods. 

3 Theoretical Background 

The methods used for extracting vector field topology are founded upon the 
theory of the qualitative behavior of dynamical systems. Most of the brief 
review of the essential theory in this section is taken from the books [DLA061 

EctHUElgHzS]. 

3.1 Dynamical Systems 

Definition 3.1 A dynamical system on E C W 1 , an open subset ofW 1 , is a 
C 1 map 4> ■ K x E — > E, where <f> = <j)(t, x) with t G M., x G E, that satisfies 

1. 0(0, x) = x \/g eE , 

2. 4>{s, x) o (f)(t, x) = cf)(s + t, x) VVeE, u ,„ eR . 

Dynamical systems are closely related to autonomous systems of differential 
equations. On the one hand, let E C W be open and let / G C 1 (i?) be 
Lipschitz continuous on E. Then the initial value problem of the autonomous 
system of differential equations 

= (3.1) 

with x(0) = xq for any xq G E has a unique solution defined for all t G K 
by virtue of the Picard-Lindelof theorem. For each initial value xq G E this 
induces a mapping (f>(t,x~o) : M x E E referred to as a trajectory of the 
system (|3.1[) . The mapping <f> then lies in C 1 (R x E) and thus is a dynamical 
system in the sense of Definition 13.11 It is called the dynamical system induced 
by the system of differential equations (|3.1[) . On the other hand, if (f>(t,x) is a 
dynamical system on E C K", then 

f(x) = -r${t,x)\ t =o 

defines a C 1 vector field / on E and for each xo G E, cf>(t,Xo) is the solution to 
the initial value problem with x = xq of (|3 . 1 [) . 
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(d) (e) (f) 

Figure 2: First-order critical points of planar vector fields classified in terms of 
the eigenvalues Ai,A 2 of the Jacobian: (a) saddle (3?(Ai )3?( A2) < 0, 9(Ai) = 
Q?(A 2 ) = 0), (b) attracting node (3?(Ai), 5ft(A 2 ) < 0, 3f(Ai) = 3(A 2 ) = 0), (c) 
attracting focus 3?(A 2 ) < 0, 3(Ai),3(A 2 ) + 0), (d) center (SR(Ai) = 

5R(A 2 ) = 0, 9f(Ai),9f(A 3 ) 7^ 0), (e) repelling node (»(Ai), 5R(A 2 ) > 0, 9f(Ai) - 
Q?(A 2 ) = 0), (f) repelling focus (»(Ai),»(A 3 ) > 0, 3(Ai),3(A 2 ) ^ 0). 

3.2 Critical Points 

An important concept in the field of dynamical systems is the notion of a critical 
point: 

Definition 3.2 An equilibrium or critical point xq £ R n of a dynamical system 
4> is a point where <j)(t,Xo) = Xq VtgR. If the dynamical system is induced by 
a system of differential equations h3.1\) . then a critical point a?o of <f> is a point 
where f(xo) = 0. If the Jacobian Df(xo) has only eigenvalues with nonvanishing 
real part, x*o is called a hyperbolic critical point. If det Df(x*o) ^ 0, then 
xq is called non-degenerate or first-order critical point. Otherwise it is called 
degenerate or higher-order critical point. 

A system of differential equations can be approximated by its linearization 
around a critical point xq without changing its qualitative behavior if xo is a 
hyperbolic critical point of that system (Hartman-Grobman theorem [Per9f j ). 
For planar linear systems, only certain types of critical points can occur and 
these can be classified in terms of the eigenvalues of the Jacobian as shown in 
Fig. [5] (here only non-degenerate cases with det Df(x~o) ^ are considered). 

3.3 Poincare Index of a Critical Point 

In order to classify critical points of vector fields one can use the notion of the 
Poincare index (or index for short) of a critical point. 

Definition 3.3 Let f = (/i,/ 2 ) T ieflC'fE) vector field on some open E CM . 
If xq is an isolated critical point of f and T C E a Jordan curve such that x~q 
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is the only critical point of f in its interior, then the Poincare index of Xq (or 
index for short) is 

I r (x-o) :=//(r) :=^J£Z, 
r 

with 9 = arctan fy- . 

ii 

It can be shown that isolated first-order critical points (i.e. isolated critical 
points for which the Jacobian of the vector field in the critical point has no 
eigenvalue of 0) have a Poincare index of ±1 and that a saddle has a Poincare 
index of —1, whereas non-saddles have a Poincare index of +1 (see Fig. [2]). 

3.4 Topological Equivalence and Sectors 

Let us now establish the fundamental notion of topological equivalence of vector 
fields: 

Definition 3.4 Suppose that f € C l {E) and g e C 1 (F) with open sets E,F C 
W l . The two autonomous systems of differential equations x — f(x) and x = 
g(x) and thier induced vector fields are said to be topologically equivalent if 
there exists an orientation preserving homeomorphism that maps trajectories of 
the first system onto trajectories of the seconds system. 

Markus |Mar54| showed that for planar C 1 systems of differential equations 
the condition of being topologically equivalent is equivalent to the systems hav- 
ing the same separatrix configurations, where a separatrix of a system (|3.ip is a 
trajectory of (|3.ip which is either a critical point, a limit cycle, or a trajectory 
lying on the boundary of a hyperbolic sector as defined below. This justifies the 
use of the term vector field topology for the topological skeleton of a vector field 
consisting of separatrices of that field. 

The notion of sectors was first introduced by Poincare [Poi93] to investi- 
gate higher-order critical points of planar systems, and later extended by Bcn- 
dixon |Ben01| and Andronov |ALGM73] , The idea is that one can describe the 
qualitative behavior of a planar C 1 vector field / in a suitable neighborhood 
N(xq) of an isolated critical point Xq of / in terms of connected regions, so 
called sectors, which form a partition of N(xq). Within each sector the trajec- 
tories of / exhibit a behavior that is characteristic for this type of sector. It can 
be shown that there exist three topologically different types of sectors: 

Definition 3.5 A sector of a critical point xq can be classified as a hyperbolic, 
parabolic, or elliptic sector according to its topological structure as shown in 
Fig.m 

Tricoche et al. [TSHOOj also used the idea of sectors to model higher-order 
critical points with a piecewise linear interpolation scheme. 

3.5 Bifurcation Theory 

Bifurcation theory is based on the notion of structural stability of a vector 
field due to Andronov and Pontryagin [AP37j: if the qualitative behavior of a 
dynamical system (|3.1|) does not change for small perturbations of the vector 
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(a) (b) (c) 

Figure 3: The three distinct topological sectors of the vector field around an 
isolated critical point with nonvanishing Jacobian (modulo the operation of re- 
versing the vector field direction): (a) hyperbolic sector (with two separatrices), 
(b) parabolic sector, (c) elliptic sector. 

field /, then that vector field is called structurally stable. If / is not structurally 
stable, the topological skeleton of the vector field changes even under small 
perturbations of the vector field / and / is said to be structurally unstable or 
to lie within the bifurcation set. 

The perturbation of the vector field is usually modeled by an additional 
parameter fj, € M: 

§ = /(^)- ( 3 -2) 

A value of \i = Ho G R for which the system p.2[) lies in the bifurcation set is 
called a bifurcation point of (|3.2p and /xo is then called bifurcation value of (|3.2p . 

Bifurcation theory has been studied extensively, see for example the book by 
Guckenheimer and Holmes [GH90] . It also explains the splitting of higher-order 
critical points into several nearby first-order critical points: it can be shown 
that, if a vector field / has an isolated critical point x*q of higher order, there 
exists a perturbation of / such that x splits into several isolated first-order 
critical points nearby. 

4 Classifying Critical Points Without Using the 
Jacobian 

In this section, we show how the Poincare index of an isolated first-order critical 
point can be computed by evaluating the sign configuration of the vector field's 
components on a finite set of sample points in a neighborhood of the critical 
point, reminiscent of the marching-cubes classification applied to isosurfaces in 
scalar fields. 

4.1 Setting 

From now on let / € C 1 (E) be a vector field / : E C IR 2 — >• R 2 defined on 
an open £cl 2 such that / is Lipschitz continuous on E and only has non- 
degenerate first-order isolated critical points. 

4.2 w-level sets, Areas of Characteristic Behavior 

We introduce the notion of areas of characteristic behavior that will enable us 
to calculate the index of a critical point. 
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Lemma 4.1 Let f = (fi, fa) be a vector field like in \4-l\ and Xq G E an 

isolated first- order critical point of f , i.e. detDf(x~o) 

=/= 0. Then, xq lies in the intersection of the 0-level sets c\ and c 2 of f\ and 
fa. Furthermore, there exists an e > such that c\ and c 2 are infinitesimally 
straight lines (i.e. lie infinitesimally close to straight lines) in an e-neighborhood 
ofx Q . 

Proof. By definition, a critical point of / has to lie in the intersection set of 
the 0-level sets of the components of /. Since / is lincarizable around xq, 
Taylor's theorem leads to f(x) ps f(x~o) + Df(x~o)(x — x~o) with x in the e-ball 
B € (a^o). Since Df(x~o) has full rank, the 0-level sets of Df(x*o) are straight lines 
intersecting in ieo. □ 

Definition 4.2 Let f = (fi,fa) T be a vector field like in \4-l\ and Xq G E an 

isolated first-order critical point of f. Then for an e > the 0-level sets c\ of 
fi and C2 of fa partition the e-ball B c (xq) of x*o into four disjoint open subsets 
A\ , . . . , Aa called areas of characteristic behavior or areas for short. 

In each area, the signs of f\ and fa do not change, i. e. for arbitrary 1 < i < 4 
and x , y G A, the following holds: 

fa / and sgn(fj(x)) = sgn(fj(y)) for j = 1,2. 

The set of areas A\, . . . , A4 around a critical point can be seen as an ordered, 
cyclic sequence of areas according to the order in which they intersect with the 
boundary of B c (xo), as illustrated in Fig. [^J We consider clockwise traversal 
direction unless otherwise noted. 

Since each of the two vector field components can either be positive or nega- 
tive in one area, there exist four different types of characteristic areas as shown 
in Fig. \5j(a). In the following, areas are written as ordered 2-tuples over the set 
{+,— } or equivalently as elements of the set {1,2,3,4}, where 1 = (+,+), 
2 = (+,—), 3 = (—,+), 4 = (—,—). Area sequences can then be written 
as ordered 4-tuples over the set of areas, i.e. as ordered ^-tuples over the set 
{1, . . . ,4}. Two areas that lie next to each other are called adjacent areas. As 
Ai and A4 are adjacent, the indices of the areas in the area sequence are cyclic, 
i.e. the area sequence A\, ... , A4 is glued together at A\ and A4. 

Remark 4.3 For two adjacent areas Ai,Ai + i of an area sequence with x G Ai, 
y G A4+1, either 

sgn(/i(x)) ^ sgn(/i(y)) A sgn(/ 2 (f)) = sgn(/ 2 (y)) or 
sgn(/i(x)) = sgn(/i(y)) A sgn(/ 2 (f)) ^ sgn(/ 2 (y)) 

holds, i.e. exactly one component flips its sign for two adjacent areas but not 
both. 

4.3 Area Sequence and Types of Critical Points 

Definition 4.4 Let xTq be an isolated first-order critical point of f defined like 
above with an area sequence 

(A±, . . . , A4). Then for a pair Ai,Ai + i of adjacent characteristic areas of the 
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Figure 4: Area sequence (Ai, . . . , A4) around an isolated critical point of a C 1 
vector field / = /2) T defined by the 0- level sets c\ and C2 of /1 and f2- The 
area sequence of a critical point a?o can be constructed by walking monotonously 
around the boundary of an e-ball B e (xo) of xq starting at an arbitrary position 
(shown as dashed line above) and collecting the intersection points of c\ and ci 
with B c (xq). 
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Figure 5: Areas of characteristic behavior and area sequences around a critical 
point C: (a) the four types of characteristic areas Ai, (b) counterclockwise and 
(c) clockwise turning behavior of the areas, yielding a critical point of index — 1 
and +1. respectively 
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Figure 6: The change of angle of the vector field for two adjacent areas; since 
only one component flips sign, the two vectors lie in adjacent quadrants I-IV 
and one has \8\ < it. 



area sequence of xq (see Fig.\5j(a)), a turning in the characteristic behavior of 
the vector field can be defined. This turning can either be a clockwise turning or 
a counterclockwise turning as defined in Figs.\^b) and\^c), respectively. Since 
Figs.\^b) and\^c) contain all 8 possible configurations of pairs of characteristic 
areas, only a clockwise or counterclockwise turning behavior can occur and the 
list in Figs.\f^b) and\Sj(c) is complete. 

Remark 4.5 One adjacent pair of characteristic areas already determines the 
whole area sequence in terms of its turning behavior. This is the case because 
the vector field components whose signs flip are alternating in the area sequence. 
Thus, an area sequence can either show a clockwise turning behavior or a coun- 
terclockwise turning behavior, but not a mixture of these and we will also speak 
of a clockwise or counterclockwise turning behavior of the area sequence as a 
whole — a clockwise or counterclockwise area sequence for short. 

We can use the turning behavior of area sequences around a critical point 
to determine its Poincare index: 

Theorem 4.6 Let f G C 1 (£') be a vector field like in \4^1\ and let xq <E E be 

an isolated critical point of f of first order, such that Df(xo) has full rank. If 
the area sequence of xq is counterclockwise, then the Poincare index of Xq is 
Ij{xq) = —1 and Xq is a topological saddle of f . If the area sequence of x*q is 
clockwise, then I ^(xq) = +1 and Xq is a non-saddle first-order critical point of 

I 

Proof. The area sequence can also be interpreted as a sequence of four qualita- 
tive samples of the vector field values lying on a piecewise linear closed curve 
r C B e (xo) that contains xq (qualitative in the sense that just the signs of the 
vector field components are sampled) . We will now prove that summing up the 
angle change of the vector field between these discrete samples and performing 
this for all four samples yields the same result as evaluating the continuous in- 
tegral of the change of angle of / over a continuous Jordan curve T around xq, 
only containing one critical point xq. Therefore, the Poincare index of xq can 
be computed by just identifying the sequence of characteristic areas around xq. 

Let / be the linearized system of / defined on an e-ball around xq. Since 
for two adjacent areas only the sign of one component of / changes, the vector 
at the first sampling point and the vector at the second sampling point lie in 
two different but adjacent quadrants around xq (see Fig- EJ) - As the value of the 
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linear approximation / changes linearly on T between two sample points and / 
is continuous, the turning behavior is uniquely determined by the two sample 
points; and since the two sample points lie in adjacent quadrants defined by the 
coordinate axes, the vector rotates about an angle 9 with \9\ < it when walking 
on r from one sample point to the next. Thus, the whole area sequence (for 
four pairs of areas) yields the total change in angle 9, with |0| < Air. As the 
area sequence wraps around the critical point, the first and the last vector of 
our sample sequence are the same and thus the change in angle of the vector 
field has to be 2kir with k e Z, see Fig.^b) and Fig.^c). As detDf(x ) + 0, 
for some e also det Df(x) ^ for x G B e (x~o), i.e. the vector field is not constant 
in a neighborhood of xq and one has k ^ 0. It is \k\ < 2 because |9| < Air. 

Thus, only k = ±1 is possible and the critical point of the linearized field / is 
of Poincare index I~(x~q) = ±1. 

Since for linear systems all critical points of index ±1 have been classified (see 
Sections 13. 21 and [3. 31 as well as Fig. [2j) , a critical point can either be a saddle (for 
Iz.{xo) = — 1), if the area sequence shows a counterclockwise turning behavior, or 

a non-saddle (i.e. a source or sink for I^xq) = —1), if the area sequence shows 
a clockwise turning behavior. According to the Hartman-Grobman theorem, 
the Poincare indices of first-order isolated critical points are invariant under 

linearization and / has a saddle in xq if and only if / has a topological saddle 
in x~q, which completes the proof of the theorem. □ 

4.4 Invariant Operations on the Area Sequence 

There are 8 possible area sequences as there are four different areas and each 
area sequence can be clockwise or counterclockwise. Since these sequences yield 
only two types of critical points distinguishable by their Poincare index, it is 
desirable to build equivalence classes of area sequences that yield critical points 
of the same Poincare index. These equivalence classes can be directly con- 
structed using elementary group theory. We refer to the book |Hup67| for a 
comprehensive overview of the theory of finite groups. 

First it is obvious but nonetheless interesting to observe that the turning 
behavior of the area sequence is invariant under rotation and also invariant under 
a simultaneous sign flip of both vector field components. Figure [5] illustrates 
this: for example, the first and the last area pair of Fig. EJb) and Fig. [5Jc) can 
be related through a sign flip, likewise the second and third area pair. 

Rotations and sign flips can be modeled as group operations. Given a group 
G and a set X, a (left) group action of G on X is a mapping o : G x X — > X 
denoted (g, x) h- > g ■ x such that e • x = x for all x £ X (here e denotes the 
neutral element in G) and (gh) ■ x = g ■ (h ■ x) for all g, h € G and all x £ X. 
The orbit of an element x 6 X under a group action of G on X is the set 
0{x) = {g ■ x : g £ G}. Orbits of a group action on a set X define equivalence 
classes on X, and the set of all orbits is a partition of X. 

The rotation of the area sequence is identical to the operation of the cyclic 
group C4 on the indices of the area sequence. The sign flip can be modeled as 
the operation of a group isomorphic to the symmetric group £2 on the signs of 
the components. From now on, the first group is referred to as the shape group 
G s and the second group as the (sign) flip group Gf. The direct product of G s 
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and Gf is called the coloring group G c — G s x Gf. 

Now we define the group action of G c on the set of area sequences. Let 
ir r : G c —> G s be the projection of G c onto G s and irf : G c — >• Gf the projection 
of G c onto Gf. Further, let g G G c be an element of the coloring group G c with 
its projections 7r r (g) = s £ G s and TTf(g) = f E G/ onto the shape group and 
the flip group, respectively. Then, g acts on an area sequence (A\, . . . , A4) as 
defined below: 

g(A 1: ...,A A ) := (/A s(1) , . . . , /A s(4) ), 

where s is a permutation of the indices {1, ... ,4} and / a self- inverse permu- 
tation on the set of areas, {1, . . . , 4}. / can be interpreted as an element that 
flips the signs of all components of the vector field: areas of type (+, +) = 1 
are mapped to areas of type (— , — ) = 4 and vice versa; an area of type 2 is 
interchanged with an area of type 3. 

The orbits of this group action on the set of all possible area sequences yield 
equivalence classes of area sequences such that each equivalence class contains 
all area sequences that can be mapped onto each other by rotations and sign 
flips. 

4.5 Interpolated Vector Fields 

Typical vector field data is given on a grid. Therefore, the vector field needs 
to be interpolated to obtain values at non-grid points. Local (i.e. cell-wise) 
interpolation schemes are often chosen for the sake of simplicity and speed. 
In the following two sections, we will have a closer look at two interpolation 
schemes: barycentric interpolation on triangles and bilinear interpolation on 
rectangles. 

Both interpolation schemes share some important properties. On the one 
hand, inside cells, both interpolation schemes are of class C°° and thus lin- 
carizable. On the other hand, they arc defined locally or in a piccewise way: 
a different interpolant is used for each cell. The interpolation is only of class 
C° across cell boundaries and it is linear and continuous along cell boundaries. 
Basic concepts and tools are developed for the simpler case of the barycentric 
interpolant. Then, these tools are adapted and extended for bilinear interpola- 
tion. 

Note that the case of linear interpolation on triangular meshes can equally 
efficiently be solved by calculating the rotation along each triangle directly. As 
there can either be no or exactly one critical point of Poincare index ±1 inside 
each cell, a cell with zero rotation has no critical point on its inside, whereas a 
nonzero rotation implies that there is a critical point inside the cell and already 
determines its Poincare index. None the less we will describe the barycentric 
setting in the following as the notions introduced there will also be used in the 
more subtle case of the bilinear interpolant. 

5 Barycentric Case 

A d-simplex s on d+1 vertices offers a natural way for the linear interpolation via 
barycentric coordinates. For the following, we will again restrict the dimension 
to d = 2, where a simplex is identical with a triangle. 
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5.1 Cells 



From now, a cell T denotes a triangle with vertices xi, xi, x 3 £ R 2 and with 
one real 2D vector (from the vector field) attached to each vertex such that a 
cell becomes a 3-tuple T = ((x 1 ,v 1 ),(x 2 ,v 2 ),(x 3 ,v 3 )) c (M 2 x R 2 ) 3 . Most of 
the time, we are only interested in the vector field values and not the position 
of the vertices, just writing T = (vi, vi, Va). The vector field inside a cell T will 
be written / : \T\ -> M 2 , when |T| C R 2 is the set of convex combinations of 
the vectors xi,X2,x 3 . Further, let / = {fi,fi) T ■ We restrict the vector field 
values on the vertices to be non-zero such that the interpolated vector field has 
isolated singularities only. 

5.2 Sector Sequences and level sets of the Interpolant 

The vector field is interpolated component-wise inside the cell using barycentric 
coordinates. The interpolation is linear, and the interpolant is continuous in 
the cell and on its boundaries. 

Definition 5.1 Given the cell T, let us look only at one component of f, say 
fi (fi)- If two adjacent vertices (i.e. two endpoints of an edge of the cell) have 
values of fi (fi) such that one value is smaller than uj and the other one bigger 
than uj, then by virtue of continuity and linearity of the interpolant on cell edges 
there exists exactly one point on the edge with an fi-value (fi-value) of uj. Such 
an edge is called an w-activc edge for f\ (fi). If the two values are both smaller 
or bigger than uj, there exists no such point and the edge is called an ^-inactive 
edge for f x (fi). 

Remark 5.2 We can now observe the following: 

1) For one cell, there exist at most two cj-active edges for each component 
and, thus, at most two points on the cell boundary with value uj for each com- 
ponent if all vertex values of each component differ from uj. 

2) Since a critical point of the interpolated vector field can only occur in a 
crossing point of the 0-level sets of the components of the vector field, a cell with 
an interior critical point requires two active edges per component. Furthermore, 
there can at most be one critical point inside a cell because the 0-level sets of 
the interpolant are straight lines. 

3) Since there can only be one critical point inside a cell, the sector sequence 
around a critical point can be found by looking for intersections of the 0-level 
sets of the components with the cell boundary. Each active edge yields exactly 
one 0-value point on a boundary edge, the position of which can be obtained 
through linear interpolation between adjacent vertices. 

If one collects the sequence of 0-value points of the components while traversing 
the boundary of the cell (as shown in Fig. [7]), one can not only construct the 
area sequence but also determine whether the 0-level sets cross and thus whether 
there is a critical point of the interpolated vector field. Let a denote a 0-value 
point of the first component and let b denote a 0-value point of the second 
component of the vector field. Then, for the sequences aabb, bbaa, abba, or 
baab, the cell does not contain a critical point, whereas for the sequences abab 
and baba, there is a critical point. These (2) = 6 intersection sequences are all 
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Figure 7: Triangle cell data structure with two edge slots on each edge. An 
edge slot is identified with a possible 0- value point of one of the vector field 
components. Sy refers to the j-th slot of edge i. 




Figure 8: Different types of critical points for the same vertex sign configuration: 

(a) saddle with area sequence (++, +— , , — +), (b) attracting node with area 

sequence (++, — h , H — ). 

possible sequences for the barycentric interpolant because there can be at most 
two 0- value points for each component on the boundary of the cell. 

The sequence of the crossings of the 0-level sets of the components together 
with the information of how the signs of the components change defines the 
sequence of characteristic areas around a critical point and thus its Poincare 
index as we have shown in Section EOl We will make use of this in the following 
section. 

5.3 An Edge-Coloring Problem 

Before we can classify critical points, a suitable data structure for describing the 
area sequence around a critical point is needed. Just looking at the sign con- 
figuration of the components in the vertices (which can be seen as a coloring of 
the vertices) does not allow us to distinguish the types of critical points because 
the area sequence is not fully determined just by the signs of the components 
in the vertices (see Fig. H]). 

We use an edge-based data structure to uniquely describe the sequence of 
0-values of the vector field components on the boundary of the cell and the 
area sequence: each cell edge is given two slots that can be filled with 0- value 
points of the components (see Fig. [7]). The slots represent the order of the two 
components' zero values while walking around the boundary of a cell. The slots 
also indicate whether a component changes from positive to negative or from 
negative to positive values as it passes through a 0- value point. There is no more 
than one 0-value for each component along a boundary edge of a cell because 
of the linearity of the interpolation along cell edges. Hence each edge can be 
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Table 1: The thirteen edge colorings. 



Edge type 


0-values of components 


& 


no 0-value, neither for fx nor for fi 


6 


(fi 


+ ->-) 


£ 3 


(h 


-->+) 


U 


{h 


+ -►-) 


& 


ih 


--►+) 


£ 6 


(fi 


+ ->- /a 


+ ->•-) 


£7 


(fi 


+ ->-, /a 


-> +) 


& 


(A 


-->+,/a 


+ ->■-) 


6 


(A 


-->+,/a 


-> +) 


60 


ih 




+ ->•-) 


£n 


{h 


+ ->-/! 


-> +) 


62 


ih 




+ ->■-) 


63 


(/a 




-> +) 



given two slots and 13 types of boundary edges can occur, as listed in Table [T] 

The notation (/1 : H — > — , fi ■ ► +) denotes that /1 changes from a positive 

value to a negative value, fi changes from a negative value to a positive value, 
and the sign change of /1 occurs before the sign change of as one traverses the 
edge on the cell boundary in clockwise direction. If a component is not listed, 
it indicates that this component has no sign change on that edge and thus no 
0- value. 

For a triangle cell, each of the three edges can now be of types £i-£i3, which 
represent a coloring of the triangle edges with 13 colors. Therefore, such an 
edge-coloring can be seen as an ordered 3-tuple over the set of the 13 colors. 
Note that not every possible 3-tuple over the set of colors is a valid coloring as 
there exist the following types of invalid colorings, i.e. colorings that are not 
(or not uniquely) mappable to a sign configuration of the vector field at the cell 
vertices: 

1) Colorings that have an invalid number of active edges: For each compo- 
nent, the cell has to have at least two active edges and the number of active 
edges has to be even. 

Example: the coloring (£2, £4, £1) would have one active edge for both f\ and fi 
and would thus be an invalid coloring. 

2) Colorings that are not mappable to a sign configuration of the vector field 
in the cell vertices at all. 

Example: a component, say fx, cannot change from + to — twice in a row. 
Thus, for example, a coloring (£2, £2) £2) would be invalid. 

Each valid coloring of a triangle can by construction be uniquely mapped to 
a sign configuration of the values of the vector field components on the vertices 
(see Fig. [TTj). However, the edge-coloring of a triangle carries more information 
than just the sign configuration of the vector field components on the vertices 
(namely a unique description of the intersection topology of the 0-level sets with 
the cell boundaries) such that several colorings of the triangle may be mapped 
to the same sign configuration. 
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Table 2: Representatives of the 8 equivalence classes of colorings of a triangle 
cell. 



Class 


Cell coloring 


Index 


1 


(1,6,9) 


-1 


2 


(1,7,8) 


+ 1 


3 


(1,10,13) 


+ 1 


4 


(1,11,12) 


-1 


5 


(2,4,9) 


-1 


6 


(2,5,8) 


+ 1 


7 


(2,11,5) 


-1 


8 


(2,13,4) 


+ 1 



5.4 Classification 

As discussed in Section l4~4l there are certain operations that leave the clockwise 
and counterclockwise turning behavior of the area sequence invariant. These 
operations can be modeled as a group action of a certain group on the set of 
area sequences. 

Since there exists a bijective map from the set of edge-colorings of a cell T 
to the set of area sequences (each coloring describes exactly one area sequence), 
one can make use of the ideas presented in Section POl to construct equivalence 
classes of cell colorings yielding the same area sequence. Then, all possible edge 
colorings of a cell can be constructed and classified by using a group action that 
builds equivalence classes of equivalent colorings. Colorings related by rotation 
or color flip thus yield the same type of critical point. 

Theorem 5.3 Every configuration of a cell T (as a coloring of the cell) that 
results in an intersection of the level sets of the two components inside ofT and 
thus contains a critical point is topologically equivalent to one of 8 representative 
configurations of which four contain a saddle first-order critical point (Poincare 
index —1) and four contain a non-saddle first-order critical point (Poincare 
index +1 ). Representatives for the 8 orbits are listed in Table® and visualized 
( Online Resource 1 ). Colorings that are not included in this list do not contain 
a critical point. 

Proof. This theorem is proven using a GAP program, see (Online Resource 3). 

The basic idea is as follows. Using a combinatoric description, all possible 
colorings of a cell arc constructed: the set of all ordered 3-tuples over the set 
of all possible colors. Then, a group action on that set is considered to build 
equivalence classes of colorings such that all pairs of elements of an equivalence 
class can be mapped onto each other by means of operations leaving the type 
of critical point invariant, as described in Section [4.41 

Let C = {£i, . . . , £13} be the set of colors as described above and let the set 
of 3-tuples over C, T = {t = (<i, . . . , t 3 ) : ti E C} be the set of colored cells or 
tiles. For the coloring group G c , which is the direct product of the shape group 
G s and the flip group G / (see Section 14. 4p , an action of G c on the set of all 
colored tiles T can be defined. Here, the shape group G s is the cyclic group C3 
as the subgroup of rotations of the symmetry group of the combinatoric triangle 
(which is given by the symmetric group S3). The color flip group is chosen to 
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be G f = ((2,3)(4,5)(6,9)(7,8)(10,13)(11,12)), which is isomorphic to S 2 ; the 
numbers in the cycles correspond to the color indices of edge colors. The choice 
of generator for Gf is obvious and unique since simultaneous flipping of signs 
of the two components results in an interchange of edges with color £2 and £3, 
£4 and £5, etc. 

Using our GAP program, first all possible colorings of a cell arc generated 
and for each orbit a representative is checked for the validity of the coloring. 
Elimination of invalidly colored cells can be done on a representative level: if 
the representative of an orbit is an invalid coloring, all other elements in the 
orbit are invalid because they are equivalent colorings. Thus, one can talk of 
invalidly colored orbits. 

After discarding all invalidly colored orbits, for the remaining orbits the 
sequence of 0- value points of the components on the cell boundary is extracted 
from the coloring. Since the 0-level sets of the components are straight lines, the 
sequence of 0-value points determine whether the 0-level sets intersect within 
the cell (thus yielding a critical point) or not (no critical point). If the 0-level 
sets intersect within the cell, the area sequence and its turning behavior are 
extracted to classify the type of critical point as saddle (Poincare index — 1) 
or non-saddle (Poincare index +1). As all critical points arc of first-order, this 
classification by the Poincare index covers all possible types of critical points. 
Since all possible colorings have been considered the list of equivalence classes is 
complete. See (Online Resource 1) for a list of all equivalence classes, including 
sample visualizations. □ 



6 Bilinear Case 

In this section, the concepts developed in the previous section are transferred to 
the bilinear case. Most elements can be immediately adopted, but some caution 
has to be taken because the interpolant is no longer linear. 

6.1 Interpolation Scheme and Cells 

Bilinear interpolation is an extension of linear interpolation for interpolating 
functions of two variables. This interpolation scheme is simple, fast to imple- 
ment, and widely used in many visualization algorithms working in a cell-wise 
fashion on uniform, rectilinear, or curvilinear 2D grids. 

A cell Q is represented by an ordered 4-tuple on the points a?i , . . . , £4 £ 
R 2 with one real-valued 2D vector (the vector field values) attached to each 
vertex such that a cell becomes a 4- tuple Q = ((afi, Vi), . . . , (2:4, #4)) C (R 2 x 
R 2 ) 4 . Mostly, we are only interested in the vector field values, just writing 
Q = (vt,.. .,Vi). 

Any uniform, rectilinear, or convex curvilinear cell can be transformed to 
the unit square [0,1] 2 in terms of a diffcomorphism, yielding local Euclidean 
coordinates s,t <E [0, 1] within the cell. Without loss of generality we therefore 
only consider the case of Euclidean coordinates in the following, i.e. the setting 
when the interpolation scheme is employed in parameter space. We use the 
following notation for bilinear interpolation: 

Let Bi : [0, 1] x [0, 1] — > R, i = 1,2, be bilinear functions for the two vector 
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components (corresponding to x and y), defined by 



Vli 




(V) 


V3i 


V4i_ 





where s, t 6 [0, 1] cl are local coordinates within one cell and v\ = («n, «i2) T ,- ■ • 
(t | 4i,t | 42) T are the values of the vector field components at the four vertices 
Then B = (Bi,B 2 ) T is called the bilinear interpolant of . . . , v 4 . 

Remark 6.1 The bilinear functions can be rewritten as 



Bi(s,t) = diS + bit + Cist + di, (6-1) 

with o, = (v 3i - vu), h = (v2i - vu), 
Ci = (vu - v 2 i - v si + vu) and d 4 = vu. 

The following properties of the interpolant will be of importance throughout 
the rest of the section: 

1. The w-level sets of B L are hyperbolas, which can be seen when writing 
the interpolants in standard form (|6.1[) and interpreting them as conic 
sections. Some examples of w-level sets of Bi are depicted in Fig. [9] 

2. The interpolation is linear in each dimension and continuous. Further- 
more, the interpolation along edges of the cell is linear: analogously to 
the barycentric case in Section [5j one can define w-active and w-inactivc 
cell edges. There exist four points on a cell boundary with value u> and 
thus most four w-active edges, if all vertex values differ from ixs. 

3. In the case = 0, Bi(s,t) is a linear function. If Cj ^ 0, Bi(s,t) is 
nonlinear and when looking at the partial derivatives J^Bi(s, t) = Cit + a; 
and j^Bi(s,t) = CjS + one sees that 2?j(s,t) has an extremal point at 
Sj(— — I — I 1 ). This extremal point is a saddle because the Jacobian DBi 
is singular at S% and for all s,t € M the Hessian Hi of Bi has a negative 
determinant (and thus its eigenvalues are of opposite sign; see Figs. |H] and 
[12] for examples). 



6.2 Bilinear Interpolation of Scalar Fields 

Let us look at the qualitative behavior of a bilinearly interpolated scalar fields 
before extending this to vector fields. We only consider one Bi here, say Bi, 
as a placeholder for a scalar field. The qualitative behavior of the interpolant 
within the cell can be described in terms of its vertex configuration: 

Theorem 6.2 Given lo £ R and a cell Q with vertex values vu, V21, U31 and 
vu, vu 7^ lo for 1 < i < 4, four cases for the qualitative behavior of lo -level 
sets of the interpolant Bi(s,t) within a cell can occur: Fig. summarizes these 
cases that depend on the classification of vertex values vu being greater or less 
than lo. Here, active and inactive edges for a cell Q are defined in the same way 
as for the barycentric interpolant in Definition \5.1\ Please note that we use the 
terminology saddle cell ( see Fig. OJ) to describe a configuration of vertex values 
Vu, which is different from the classification of a critical point as a saddle point. 
In the following, we denote the first always by the compound term saddle cell. 
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(i) («) I 



A 
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(iii) 
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(iv; 
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Figure 9: Examples of the four qualitatively distinct configurations of 0- level 
sets of bilinearly interpolated scalar fields: (i) inactive cell, (ii) single active cell, 
(iii) couple active cell, (iv) saddle cell. Bounding boxes (see Lemma HT4| arc 
shown in gray. 



Table 3: Different cell types by vertex sign configuration. 



Cell type 


Vertex sign configuration 


inactive 


(----),(+,+,+,+) 


single active 


(---+),(+,---),(+,+,+,-), 
(-+,--),(--+,-),(-+,+,+), 
(+,-,+,+),(+,+,-,+) 


double active 


(-,-,+,+),(+,-, -,+),(+,+,-,-), 
(-+,+,-) 


saddle cell 


(-,+,-,+),(+,-+,-) 



All of the statements follow immediately from the continuity of the inter- 
polant and its linearity along cell boundaries. 

Remark 6.3 For an easier notation of the different cell types in terms of the 
definitions of Theorem 16.21 we write a cell Q = (#1,^2, 753, #4) as a 4-tuple over 
the set {+, — } where the i-th entry of that tuple is set to + if vi > lj and to 
— otherwise. Table [3] summarizes this tuple notation, referred to as vertex sign 
configuration, as in the barycentric case. 

Lemma 6.4 (Bounding-box lemma) Let Bi(s,t) be a bilinear interpolant 
on a rectangular grid cell Q and let the cell be 0-active. For each curve of the 
0-level set connecting two 0-value points A(A x \A y ) , B(B x \B y ) on the boundary 
of the cell, the Cartesian product 

L = [mia{A x ,B x },max{A x ,B x }] 
x [minjAj,, By}, max{A y , B y }] 
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defines a bounding box: the curve of the 0-level set of Bi is inside L except at 
the O-value points at the border of the cell. See Fig. \Qfor an illustration of the 
set L for some sample configurations. 

Proof. We only give a sketch of the proof. The basic idea is the following. 
Let the bounding box L be given as the convex hull of the four points 

Pi (mm{A x ,B X }\ min{ A y ,B y }), 
P 2 (min{ A,, ,B X }\ max{ A y ,B y }), 
P 3 (m&x{A x , B x }\ max{A y , B y }), 
P 4 {m&x{A x , B x }\ min{ A y , B y }). 

Note that A,BeL. Then, the bilinear interpolant Bi(s,i) restricted to the 
bounding box L can be expressed in terms of another bilinear interpolant Bi(s, t) 
defined on L and interpolating the values of the points Pi, . . . , P4. To leave L, 
a 0-isoline of B\ has to cross the boundary of L. Thus, a crossing 0-isoline can 
only occur where the value of Bi on the boundary of L is 0. One can then show 
that apart from the two 0- values in the corner points A and B of a L, there 
exist no other boundary points of L with 0-valuc. Therefore, such a crossing 
cannot occur and the 0-level set of Bi completely lies inside L and cannot leave 
L. 

Since the interpolation function Bi along the boundary of L is continuous 
and linear, it is sufficient to show that only two points in the set L have 0-value. 
It then follows immediately that these are the only points on the boundary with 
zero value. This can be shown easily for a representative of each of the cases 
listed in Theorem IO and Fig. [SJ □ 

6.3 Analytical Computation of the Position of Critical 
Points 

The position of critical points of the bilinear interpolant B = (Bi,B2) can be 
computed as intersection points of 0-level sets of B\ and B2 by solving 

ui +\ ( \ ( a x s + bit + c x st + di \ * 

B M=\ B 2 (s,t) ) = { a 2 s + b 2 t + c 2 st + d 2 ) =°- ^ 

Since both B\(s, t) = and B 2 (s, t) = are linear or quadratic equations in s, t, 
the system can either be combined into a single linear or quadratic equation in s 
or t. We will refer to this equation as the intersection equation. Since the linear 
case occurring here can be treated in exactly the same way as the barycentric 
case described in Section^ this degenerate case will be excluded in the following 
and we will only speak of the (quadratic) intersection equation with discriminant 
A. The equation can cither have no real solution (A < 0), two different real 
solutions (A > 0), or a double real solution (A = 0), and the number of solutions 
of the intersection equation gives the number of critical points, see Fig. 1101 

Since the sign of the discriminant A determines whether the vector field has 
no, one, or two critical points, it can be interpreted as a bifurcation parameter 
of a dynamical system, where A = is the bifurcation value. Also A can be 
interpreted as a measure for the stability of the critical point — critical points 
with A = are generally not stable, whereas the case with two (A > 0) or 
no intersection points (A < 0) can be more or less stable according to the 
magnitude of |A|. 
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Given a rectangular cell with the vector field values 
V\ , . . . , V4, one can write the components of the bilinear interpolant in normal 
form like done in (|6.2[) (see Remark |6. II) . 

The intersection points of the 0-level sets of the components are obtained 
by solving (|6.2[) . First, we determine which of the interpolants is linear in s 
and t (i.e. for which a = 0); if both are linear, (|6.2[) becomes a linear system 
and can be solved separately yielding one or no solution (and not infinitely 
many solutions as the singularities arc restricted to be isolated). If one of the 
interpolants is linear and the other nonlinear, the linear one can be solved for s 
or t and this can be plugged into the other to yield a quadratic equation in s or 
t. If both interpolants are nonlinear, one can be solved for s or t and plugged 
into the other. 

All these equations have at most two real solutions yielding one coordinate of 
the 0-value points of the interpolated vector field. The other coordinate can be 
obtained through the linear or quadratic equation between s and t from the first 
step. A solution within the cell is found by restricting the coordinates to be in 
[0, l] 2 . Note that the vector field values for each component are always restricted 
to be nonzero at cell vertices and that the two components are assumed not to 
be exactly the same. 

6.4 Level Sets and Sector Sequences 

Each solution of (|6.2[) lies in the intersection set of the 0-level sets of B\ and B2 ■ 
The following cases can occur when looking at the intersections of level sets of 
B\ and B2: 

Theorem 6.5 Let c±, C2 be the 0-level sets of B\, B2 and let A e M be the 

discriminant of the intersection equation for a cell Q. Then, the following cases 
can occur: 

1. A < C\, C2 are disjoint and there is no critical point (see Fia MOV a)), 

2. A = •<=>■ c\, C2 touch in one point (see Fig. \10Y b) ), 

3. A > •<=>• C\, C2 intersect in two points (see Fig. \10V c)). 

Proof. Clearly one has: c\ and C2 are disjoint A < which proves 1). The 
same holds for 2). Now it remains to show 2), i.e. that c\ and C2 cannot touch 
twice and that the intersection equation has a double solution iff c\ and C2 touch 
in one point. 

Suppose that c\ and C2 touch in one point as shown in Fig. HOf b) . A small 
perturbation of the vector field values now either results in a case where the 
hyperbola branches do not intersect (i.e. A < 0) as shown in Fig. ITUf a). or it 
yields a case where c\ and C2 intersect in two points, i.e. where the intersection 
equation has two different solutions as shown in Fig. UOf c). As each pair of 
branches can be perturbed separately a double touching case cannot exist as 
this could be perturbed to yield a triple intersection of the branches which is 
a contradiction to the maximum number of real solutions of the intersection 
equation. This proves 3). □ 
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Figure 10: Intersection cases for the 0-level sets (red and green lines) of the 
two vector field components: (a) no intersection with A < 0, (b) one touching 
point C with A = 0, and (c) two intersection points C\ and C2 with A > 0. 
The saddle points Si (i = 1,2) of the two interpolants Bi and their asymptotes 
(dashed lines) are depicted in red for B\ and green for B 2 - 

6.5 Critical Points and Area Sequences 

Let us now look at the types of critical points yielded for the different intersection 
cases for the 0-level sets of the components as investigated in Theorem IG.5I 

1) The case A < 0, where the 0-level sets do not intersect, is trivial. 

2) For the case A > 0, when the 0-level sets of the components intersect 
twice, yielding two first-order isolated critical points of the interpolated vector 
field of first order as none of the directional derivatives is in the intersection 
points. For each critical point separately, the same things hold as for critical 
points in the barycentric case; in particular we can look at area sequences, 
turning behavior of those, etc. for each critical point separately. 

Remark 6.6 The area sequences of the two critical points always share three 
characteristic areas and these are traversed in opposite directions for each of 
the points (i.e. clockwise turning for an area pair becomes a counterclockwise 
turning and vice versa) . Therefore, one of the two points has a Poincare index 
of —1 and the other one has a Poincare index of +1. Another way to see this 
is the fact that to first observe that there are at most two critical points inside 
a cell and that the Poincare index of a cell can either be —1, or +1. Now, 
by virtue of the summation theorem for the Poincare index, two critical points 
inside a cell of the same index would force the cell to have an index of ±2, a 
contradiction. 

3) For the touching case (A = 0) of the 0-level sets, the tools developed in 
the previous sections cannot be applied: when linearizing the field around the 
touching point, the 0-level sets of the two components are tangential to each 
other and the Jacobian has at least one zero eigenvalue. Thus, the critical point 
is not of first-order and a well-defined area sequence does not exist for this case. 
However, a perturbation of the field yields two critical points, one of index +1 
and one of index — 1. Then, Poincare index theory states that the one critical 
point yielded by the touching of the 0-level sets has an index equal to the sum 
of the indices of the two critical points it splits into, i.e. the second-order critical 
point from a touching of the 0-level sets has index 0. 
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6.6 Edge Coloring and Invariant Operations on the Col- 
orings 

As the bilinear interpolant is linear along the cell edges, one can employ the 
same ideas as in the barycentric case in order to determine existence and types 
of critical points inside a cell, again using an edge-based "slot" data structure. 
Using a data structure with two slots for each edge of a cell, one can see this 
as a cell coloring problem with a coloring of a cell represented as a 4-tuple over 
the set of 13 colors like described in Section [531 

Investigating the intersection topology of the 0-level sets becomes more com- 
plicated because the 0-level sets are no longer straight lines. Using the same 
notation as for the barycentric case, a sequence of 0-value points aabb (where 
again a stands for a 0-value point of the first component of the vector field and 
b for a 0-value of the second component) not necessarily has no intersections. 
Here, it can either yield a case where the isolines have no point in common, 
touch in one point, or intersect twice as shown in Fig. [TJ 

The question is whether for a sequence of 0-value points of the form . . .aa . . . 
or . . . bb . . . the subsequence aa or bb can yield a touching or intersection of the 
0-level sets or not. If not, the part aa or bb of the sequence is irrelevant for 
the intersection topology and can be "collapsed", e.g. a sequence baababab for 
which aa cannot yield a touching of the isolines could be reduced to bbabab. If 
now bb would yield no intersection we could reduce once more to yield abab, a 
sequence which cannot be collapsed any more. 

For a fully reduced sequence, we can either determine the number of cross- 
ings and thus critical points inside the cell as done in the barycentric case (for 
sequences of the form (ab) n , n € N ) or we know that the case is value- depen- 
dent, i.e. the number of intersections of the isolines cannot be determined just 
by looking at the topology of the 0-value points on the boundary of the cell. 
The latter is the case when the fully reduced sequence contains a subsequence 
aa or bb. 

The question remains: how can one determine which subsequences of type aa 
or bb can be collapsed? Let us first focus on non-saddle cells. The bounding box 
lemma fLemma l6.4[) gives us valuable information for the case of single or double 
active cells. It states that cases aabb with single active cells for both components 
lying at opposite vertices of one edge cannot intersect, whereas all other config- 
urations with single and double active cells for the two components — although 
having no intersection in the topology of 0- value points on the boundary — can 
be perturbed to yield a touching or an intersection of the isolines and thus are 
value-dependent. 

6.7 Saddle Cells 

There are several subtleties that have to be considered when dealing with saddle 
cells. For one or both components of the vector field, both branches of the 0- 
level sets intersect with the cell. The topology of the 0-level sets inside the cell is 
not uniquely determined by the sequence of 0-value points of the components on 
the cell boundary. This ambiguity can be resolved by adopting the asymptotic 
decider }NH91j . 

To determine the intersection topology for the 0-levelsets inside the cell, one 
can sort the 0-value points by their local coordinates inside the cell in ascending 
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Figure 11: Hierarchy of data structures to describe the sequence of 0- value 
points on the cell boundary from fine to coarse: (a) double-edge data structure 
consisting of two double edges with four slots, which can be projected via tt\ to 
(b) the edge data structure consisting of four edges with two slots each that can 
finally be projected via -ki to (c) the vertex sign configuration of the components. 



(or descending) order. Then, for each component, the first and second and the 
third and fourth point in the sorted sequence of 0-value points belong pairwisc 
to the same branch of the 0-level set. 

To model this behavior, the previously used edge-based data structure is 
insufficient because the asymptotes imply a "linking" of two opposite edges 
of a cell. Therefore, the data structure needs to be extended to describe the 
position of the 0-value points in the sequence sorted by local coordinates. A 
natural choice for this is a double-edge data structure shown in Fig. lllf a). which 
describes a cell by a pair of double edges with four slots for each edge. When 
filling the slots with 0-value points of the components, their position in the 
sorted sequence for each coordinate is the index in the respective double edge if 
one arranges the slots as in Fig. [TlT aV 

Each valid double-edge configuration can be uniquely mapped to a valid 
edge coloring configuration via a natural projection 7Ti (every valid double-edge 
configuration can be uniquely mapped to a edge coloring of a cell), which in turn 
can be mapped to a vertex sign configuration via a natural projection ~ki (every 
valid edge coloring can be uniquely mapped to a vertex sign configuration of a 
cell), as described before for the barycentric case. Figure fTTI illustrates all three 
levels of representation. 

Using the double-edge data structure, one can now construct all possible 
double-edge configurations for the case of a saddle cell and check for the in- 
tersection topology of the 0-level sets of the interpolant. As we will see in the 
following, the asymptotic decider plays no role in the types of critical points 
yielded inside a cell, i.e. there is no need to pay attention to how the 0-value 
point sequence is reduced when determining the intersection topology and area 
sequence for saddle cells. 

An open question is how the 0-value sequence can be reduced in a valid 
way as it was done for non-saddle cells before. The tools provided by the 
bounding box lemma (Lemma 16.4(1 and Theorem 16.51 are instrumental here. 
Figure [12] shows models of the two cases that contain reducible subsequences of 
type aa or bb in their sequences of 0-value points on the boundary. All other 
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Figure 12: Examples of topologies of 0-value points of the components on the 
boundaries that can be reduced: (a) baababba, which can be reduced by col- 
lapsing the branch in the top right which cannot touch or intersect with other 
branches (bounding box lemma) to bbabba, (b) abababba, which can be reduced 
to abab as the two branches in the top left cannot touch or intersect with others 
(bounding box lemma) or themselves ( Theorem I6.5|l . Again, the saddle points 
Si (i = 1, 2) of the two interpolants Bi and their asymptotes (dashed lines) are 
depicted in red for B\ and green for B 2 - 




Figure 13: Two example visualizations of vector fields for which the 0-value 
sequence on the cell boundary cannot be reduced, yielding two critical points 
inside the cell. 



configurations can yield a touching or double intersection of the branches for the 
same boundary topology and thus are value- dependent (by virtue of Lemma 16.41 
and Theorem I6.5|) . 

Now that we have the tools to construct and reduce 0-value point sequences 
for cells with bilinear interpolants, let us proceed in the same way as we did for 
the barycentric case, classifying all cell colorings in terms of the types of critical 
points they (may) yield inside the cell. 

6.8 Classification 

Theorem 6.7 Every configuration of a cell Q (as a coloring of the cell) that 
yields an intersection or touching of the level sets of the two components inside 
of Q is topologically equivalent to one of 74 representative configurations, see 
( Online Resource 2) for a complete list: 
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• 38 configurations have exactly one first-order critical point, the index of 
which is ±1 and which is determined by the sequence of characteristic 
areas on the boundary of the cell. 

• 4 configurations have two critical points, one of which is a saddle (i.e. 
index —1) and the other a non- saddle (i.e. index +1). Both critical points 
are stable. See Figure \TB for visualizations of two sample cases. 

• 32 configurations are value-dependent, i.e. they have either no, one, or two 
critical points. The case of a single critical point is an unstable bifurcation 
point of the underlying dynamical system. 

Colorings that are not included in this list cannot have a critical point inside 
the cell. 

Proof. The proof is analogous to the proof of Theorem 15.31 The corresponding 
GAP program performs the following steps. 

1) Create all possible cell colorings with the 13 colors and use a group ac- 
tion of the coloring g roup Gc — ■ G $ x G f given by the direct product of the 
shape group G s and the flip group Gf on the set of all colored cells to obtain 
equivalence classes of colored cells (that contain the same number and types of 
critical points). The flip group only depends on the colors and as these have not 
changed it can be chosen as in the proof of Theorem 15.31 The shape group G s 
is chosen as the cyclic group C4 of the rotations of a cell by -| as a subgroup of 
the full symmetry group of the square, which is the dihedral group D 8 of order 
8. 

2) Examine the representatives of the orbits: sort out invalid colorings as 
in the proof of Theorem 15.31 and determine the intersection topology and the 
reduced 0-value sequence of the 0-level sets of the two components using the 
properties described in the preceding sections. This yields an area sequence 
that can be used to determine the number and types of critical points inside the 
cell. Here, four classes of configurations are distinguished: cases with exactly 
one critical point, cases with exactly two critical points, cases that are value- 
dependent, and cases that yield saddle cells for any component. 

3) Determine the turning behavior of the reduced area sequence and thus 
the type of critical point for configurations with exactly one critical point. Va- 
lue-dependent and double intersection cases need no further treatment as the 
number and types of the critical points are not fully determined by the area 
sequence. 

4) For saddle cells: construct all possible double-edge configurations and the 
reduced 0-value sequence for each configuration yielding an area sequence. 

As it turns out, for saddle cells all double-edge configurations yield the same 
number and the same types of critical points such that the asymptotic decider 
indeed does not influence the types of critical points inside a cell as claimed in 
Section I5T71 For the other cases, the number of equivalence classes given in the 
theorem arise. A full list of cases, including example diagrams, is provided in 
(Online Resource 2). □ 
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Figure 14: Forming a super cell consisting of the quadrangular cells Qi and Q2 
for a case where a critical point C lies on the cell boundary of Q\ and Q2- The 
super cell is then triangulated by T\ , . . . , T 6 and the vector field (shown in gray) 
is interpolated barycentrically inside each triangle Tj. 



7 Boundary Points and Closed Streamlines 

So far the issue of critical points on cell boundaries has not been addressed and 
is somewhat delicate because the interpolation scheme is only C° across cell 
boundaries for both the barycentric and the bilinear cases. 

To resolve this issue, we employ a cell clustering as described by Tricoche 
et al. |TSH00j : the current cell with a critical point at the cell boundary and 
its neighboring cell sharing the edge with the critical point are clustered into 
a super cell, see Fig. [TJJ This process is iterated as long as there are critical 
points on the cell (or super cell) boundary. Then, a new (artificial) critical 
point C is positioned inside the super cell at a position given by the averaged 
positions of the critical points inside the super cell. This new critical point C 
can be of arbitrary order and complexity. In order to determine the topology of 
the vector field inside the super cell, the super cell is triangulated with a new 
inner vertex at the position of C . Subsequently, a piecewise linear vector field is 
constructed inside the super cell as the union of the barycentrically interpolated 
vector fields on the triangles, using the information from the boundary of the 
super cell. As the direction of the new vector field does not change on rays 
emerging from C , hyperbolic sectors around C can be identified by looking at 
the configuration of the vector field on the boundary of the super cell. After 
identification of the boundaries of hyperbolic sectors, the separatrices can be 
traced beginning from the super cell boundary in the usual way. Figure [T5l shows 
how the different sector types around the critical point C can be determined 
by using a sequence of orientations of the vector field in relation to the position 
vector on the boundary of the super cell with respect to the origin C . We refer 
to [TSH OOl for details of the method. This is a local method that does not alter 
the vector field behavior outside the super cell; thus, it is a robust and simple 
way to deal with critical points on cell boundaries. Note that the data of the 
randomly generated test vector fields as shown in Table |4] suggest that boundary 
critical points occur extremely seldom such that the additional steps needed by 
the algorithm in that case are negligible regarding its speed. 

Another problem not addressed so far is the task of finding limit cycles of vec- 
tor fields. Since this problem has been addressed in various publications already 
(see |WS01l ITWHS04| ) and the methods presented there can be incorporated 
in our algorithm, we will not deal with this case here. 
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Figure 15: Possible configurations of the orientation of the vector field in rela- 
tion to the position vector with respect to the origin C. P denotes a parallel 
configuration, O an orthogonal one. The appended signs + and — distinguish 
the two possible directions of the vector field w.r.t. the position vector in each 
case. 

8 Algorithm 

The algorithm consists of three main stages. First, all cells that may contain 
critical points are identified. Then, these cells are examined more thoroughly 
to check whether they contain one or more critical points and whether cell 
clustering is needed. This step is performed by computing the coloring of the 
cell and using a lookup table to determine the class of the cell coloring as 
defined in Theorems 15.31 and 16.71 Finally, the positions of the critical points 
inside the cells are computed analytically. For non-saddles, the type can either 
be classified using the Jacobian at that point or the sector-based idea described 
in |TSH00| . which is also used to determine the types of higher-order critical 
points. Note that the exact type of first-order non-saddles is not required to 
build the topological skeleton, which is only based on trajectories originating 
from hyperbolic sectors. 

8.1 Generic Algorithmic Skeleton 

1. Traverse cells, mark cells that have at least two active edges for each 
component by looking at the sign configuration of the vector field at the 
vertices. 

2. Compute cell colorings of the marked cells and fetch the configuration 
class for each cell in the lookup table; a perfect hash function can be used 
for this. Deal with critical points on the boundary of the cell using the 
cell-clustering approach as described in Section [JJ 

3. Compute the location(s) of the critical point (s) inside the cell. For the 
barycentric linear system is solved. For the bilinear case, one 
can write the intcrpolant for each component in normal form Bi(s,t) = 
ttjS + bit + Ci-st + di and combine the two equations into the linear or 
quadratic intersection equation in s or t. Then the discriminant A of 
the intersection equation and its solutions are computed. For both cases, 
solutions lying outside the cell are dropped. 
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4. Construct the topological skeleton of the vector field in the usual way 
tracing streamlines at the boundary of hyperbolic sectors and connecting 
these to other critical points or boundary points when they get close. 
The boundaries of hyperbolic sectors can either be found as described 
by Tricoche et al. |TSH00| or by using the eigenvectors of the Jacobian 
at the critical point. Our implementation employs an adaptive explicit 
Runge-Kutta method for streamline tracing. 

8.2 Performance and Topology Simplification 

We propose that, in order to overcome numeric problems in the vicinity of 
critical points, streamlines are only traced from and to boundaries of cells that 
contain critical points. Within these cells, streamlines are approximated by 
a straight line connecting the intersection point of the trajectory and the cell 
boundary with the critical point. This approach has the advantage that typically 
many small steps of an explicit solver in the vicinity of the critical point can be 
saved, thus speeding up the calculation of the topological skeleton. 

Further simplifications leading to speed-ups can be considered. First, for 
small cells, the location of the critical point can be reasonably well approximated 
by placing the critical point at an arbitrary position — say, the center — of the 
cell. Second, the topology of the vector field for cells with two critical points can 
be simplified by replacing the two first-order critical points by a single artificial 
second-order critical point. 

The position of the second-order critical point can be approximated by the 
midpoint of the straight line connecting the two first-order critical points. 

8.3 Numerical Stability 

To be numerically stable an algorithm has to yield consistent results for the 
same input data, no matter how the input data is ordered. This is especially 
important for floating-point data, where the result of interpolation may depend 
on the interpolation direction and order of instructions. For the algorithm in 
Section 18.11 this means that no matter how a cell is oriented in terms of the 
geometric location of its vertices, the result of the calculations has to be invari- 
ant for the two cells. This can be achieved by choosing a unique interpolation 
direction along cell edges such that this choice always yields the same sequence 
of calculations on the floating-point vector field data. Then, by virtue of con- 
sistency of IEEE floating-point arithmetic |IA85j . the interpolation result and 
the corresponding vertex and edge classifications have to be identical. 

Consider a cell with vertices Pi, ... , P n , Pi{xi\yi) 
(i = 1 . . . n), where Xj and yi are floating-point numbers. 

Then, without loss of generality one can choose the interpolation direction along 
an edge between two vertices to always point from the vertex with smaller y- 
coordinate to the vertex with greater y-coordinate. If the two y-coordinatcs 
agree, then the same argument can be applied to the x-coordinates of the ver- 
tices. This method always yields a unique interpolation direction as: (1) the 
inequality comparison operator for IEEE floating-point numbers is consistent in 
the sense that if fi and f% are two floating-point numbers, f\ > f% => fi < /i; 
(2) the equality comparison operator for IEEE floating-point numbers is com- 
mutative, i.e. 
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h = h => h = .fi- 
Our algorithm and the classification of the critical points depend on a con- 
sistent orientation of traversing the boundary of a cell. This orientation can be 
defined by the sign of the oriented area of a cell. For a triangle with vertices 
P\{x\\y\), P<i{xi\y?), P 3 (x 3 \y 3 ), the oriented area is 

A* = -{x 2 - £1X2/3 - Vi) - (V2 - yi)(x 3 - x-i). 

For the case of quadrilateral cells, a cell is split into two triangles and the 
oriented area sequence of one of the two triangles can be used. 



9 Results 




Figure 16: Vector field topology of the oceanic flow data set. Black dots indi- 
cate detected critical points, blue lines show the topological skeleton, annota- 
tions provide the classification (AN: attracting node, RF: repelling focus, SAD: 
saddle) . (a) Incorrect topology obtained by a method that intersects linearized 
0- level sets of the x-component (red) and the y-component (green). Notice how 
trajectories terminate in critical points that are not detected as such in the top 
left of the image, (b) Correct topology obtained with our method. The pair of 
critical points inside a value- dependent cell is correctly identified (framed area 
shown magnified in the bottom right). 



We illustrate the results of our algorithm for four test cases. The first test 
case consists of randomly generated bilincarly interpolated vector fields: for 
each cell vertex, the two vector field components are random numbers lying 
in [-1,1]\{0}, quantized to 10" 6 . Machine accuracy is set to e = 10" 9 and the 
threshold of a critical point being of second-order is chosen to be A < 0.05, where 
A is the discriminant of the intersection equation as described in Scction l5751 We 
have generated and examined 10 s cases. The second test case is a CFD dataset 
of an oceanic flow in the Baltic Sea (simulation courtesy of Kurt Frischmuth, 
University of Rostock) , see Figure [T()J The data set is given on a uniform grid 
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Figure 17: Vector field topology of the air flow data set calculated with our 
algorithm. Scparatrices are shown in blue along with a glyph-based visualization 
of the vector field, where the arrows encode the vector direction and the vector 
norm is color-coded. 

of size 100 x 112, and the vector field is rotated by 90 degrees, c.f. Theisel et 
al. |TRS03] . 

The third data set examined with the algorithm is a planar slice of a simu- 
lated unsteady flow of air in a hot room at a fixed time step of the simulation 
(data courtesy of Filip Sadlo, Universitat Stuttgart). See Figure [T7] for an illus- 
tration of the flow topology extracted with our algorithm. Again, the data set 
is given on a uniform grid, this time of size 100 x 100. 

As the visualization of flows on surfaces has recently become a focus of 
research in the field of topology- based vector field visualization [LCJK+09] . 
we chose as a fourth test case a tangential vector field on a torus given on a 
uniform grid of size 50 x 50 in parameter space and computed its topology with 
our algorithm. 
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Figure 18: Example of a vector field on a torus. Left: torus immersed in 
Euclidean 3-spacc, right: torus in parameter space. The top and the bottom 
as well as the left and the right border are identified, respectively. Again, the 
vector field is also visualized via glyphs. 
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Table 4: Occurrences of critical points in test data sets. Relative occurrences are rounded to three decimal places. 



Class 


Random field 


Ocean flow 


Air flow 


Torus 




absolute 


relative 


absolute 


relative 


absolute 


relative 


absolute 


relative 


Overall 


Topology-dependent 


62887918 


0.629 


11089 


0.990 


9972 


0.997 


2492 


0.997 


Value-dependent 


37112082 


0.371 


111 


0.010 


28 


0.003 


8 


0.003 


Boundary point 





0.000 





0.000 





0.000 





0.000 


Number of critical points 


None 


64094290 


0.641 


11141 


0.995 


9984 


0.998 


2496 


0.998 


At least one 


35905710 


0.359 


59 


0.005 


16 


0.002 


4 


0.002 


Exactly one 


33202155 


0.332 


57 


0.005 


16 


0.002 


4 


0.002 


Exactly two 


2087768 


0.021 


2 


0.000 





0.000 





0.000 


Topology-dependent cases 


No critical point 


29685763 


0.472 


11032 


0.985 


9956 


0.996 


2488 


0.995 


Saddle 


16600916 


0.264 


33 


0.003 


4 


0.000 


2 


0.001 


Non-saddle 


16601239 


0.264 


24 


0.002 


12 


0.001 


2 


0.001 


Saddle & non-saddle 





0.000 





0.000 





0.000 






Value- dependent cases 


No critical point 


34408527 


0.927 


109 


0.010 


28 


0.003 


8 


0.003 


Saddle & non-saddle 


2087768 


0.056 


2 


0.000 





0.000 





0.000 


Higher-order 


615787 


0.017 





0.000 





0.000 





0.000 



Table S] documents the statistics with regard to number and types of critical 
points for all test cases. One interesting observation is that, although higher- 
order critical points might occur, their frequency is very low (of course depending 
on the threshold value for A). On the other hand, two critical points within a 
cell are (at least in the case of the oceanic flow) more common and need to be 
considered in practice. The most interesting finding is that the realistic ocean 
and air flow data sets contain only very few value-dependent cases (less than 1 
percent). Therefore, our algorithm can in both cases identify and classify more 
than 99 percent of the critical points by just analyzing the topology of 0- value 
points on the cell boundary — without the need to evaluate the Jacobian, solve 
a quadratic system, or apply the sector method |TSH00j . Critical points on the 
boundary did not occur in this data series such that one can assume that these 
occur extremely seldom. Even if the quantization of the values for the random 
vector field test case is decreased to 10 -4 , only 46 of the 10 s generated cases, 
i.e. 0.000046%, had critical points on cell boundaries. Generally speaking: The 
coarser the quantization or numeric accuracy of the computations the likelier 
the occurrence of critical points on cell boundaries becomes. None the less we 
believe their occurrence to be very seldom in practice. 

Figure [in] shows qualitative results for the ocean data set, illustrating the 
difference between a linearized version of the vector field and the original bilinear 
version: Fig. ITBT a') shows the (incorrect) topological skeleton obtained from the 
linearized vector field and Fig. I16f b) shows the correct version produced by our 
method. The differences arise for value- dependent cells that contain two critical 
points missed when linearizing the 0-isolines of the vector field's components. 

In terms of flows on surfaces an example in form of a vector field on the 
torus is examined. The standard 2-torus T 2 can be parametrized via 

/(u, v) = ((a + fecosit) cosw, (a + 6 cos it) sin?;, bsinu) , 

with u, v, a, b € R, < u, v < 2ir, < b < a. Using this paramctrization and its 
Jacobian Df one can project two-dimensional vector fields defined in parameter 
space [0, 2ir) 2 C M 2 and their topology obtained with our algorithm to the torus, 
sec Figure [TBI 

10 Conclusions and Future Work 

We have presented a novel approach to finding and classifying critical points 
according to their Poincare index for barycentrically and bilinearly interpolated 
vector fields on triangular and rectilinear grids in parameter space, respectively 
Our approach is cell-based and efficient through the use of lookup tables remi- 
niscent of the classification of scalar fields by the marching cubes algorithm. Our 
algorithm is able to deal with critical points on cell boundaries, to detect second- 
order critical points for bilinearly interpolated vector fields (together with pro- 
viding a measure for the stability of such critical points) , and to simplify vector 
field topology inside a cell by substituting two first-order critical points with 
one of second order. We have demonstrated that, for practical applications, the 
number of critical points and their Poincare indices can be identified by just 
examining the intersection topology of the 0-level sets of the interpolants with 
the cell boundaries. We have put our algorithm on a sound mathematical foun- 
dation and showed how group theoretic tools and a combinatoric description via 
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cell colorings can be used to solve the problem of critical-point classification, 
which may not seem to be of a combinatoric nature at first glance. 

As already pointed out in Section [31 the topology-based visualization of flow 
on surfaces has become a prominent field of research during the last years (c.f. 
LCJK+09] ) and it thus would be of interest to be able to apply the algorithm 
presented in this paper to vector fields on triangulations and quadrangulations 
of arbitrary surfaces. A simple proof of concept was given in the case of the 
torus in Section |H1 Note that in the quadrangular case the algorithm is well 
suited to be used in combination with the surface paramctrization algorithm 
QuadCovcr by Kalberer et al. |KNP07| due to its structure and it thus seems a 
natural step to combine these two methods — this is research in progress. 

More fundamental open questions for future work include: can this method 
be extended to higher dimensional domains or to other interpolants such as 
tensor-product cubic? More generically, it might be of interest for other areas of 
visualization to see if they might benefit from combinatoric and group-theoretic 
approaches. 
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